Get cells of Belgian EEA reference grid containing protected areas

Goal

Given the EEA reference grid of Belgium at 1km resolution and the Belgian Natura2000 protected areas, this document assesses the cells covering the protected areas.

Setup

library(tidyverse)  # To do datascience
library(tidylog) # To provide feedback on dplyr functions
library(here) # To find files
library(sf) # To work with geospatial data
library(BelgiumMaps.StatBel)  # To load Belgian administrative boundaries
library(INBOtheme)  # To apply INBO style to plots
library(leaflet)  # To make interactive maps

Read input data

Import EEA reference grid

We import UTM grid data of Belgium at 1 by 1 km resolution. All grids have CRS (Coordinate Reference System) equal to Belge Lambert 1972:

be_grid <- st_read(here::here("data", "external", "utm1_bel"))

## Reading layer `be_1km' from data source `C:\Users\damiano_oldoni\Documents\INBO\repositories\pipeline\data\external\utm1_bel' using driver `ESRI Shapefile'
## Simple feature collection with 51726 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 3768000 ymin: 2926000 xmax: 4080000 ymax: 3237000
## CRS:            3035

be_grid_crs <- st_crs(be_grid)
be_grid_crs

## Coordinate Reference System:
##   User input: 3035 
##   wkt:
## PROJCS["ETRS89 / ETRS-LAEA",
##     GEOGCS["ETRS89",
##         DATUM["European_Terrestrial_Reference_System_1989",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6258"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.01745329251994328,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4258"]],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     PROJECTION["Lambert_Azimuthal_Equal_Area"],
##     PARAMETER["latitude_of_center",52],
##     PARAMETER["longitude_of_center",10],
##     PARAMETER["false_easting",4321000],
##     PARAMETER["false_northing",3210000],
##     AUTHORITY["EPSG","3035"],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH]]

Read Belgian protected areas Natura2000

GIS data of protected areas as defined in Natura2000 can be found on related webpage of European Environment Agency. We downloaded the GIS data as shapefile in folder data/external/protected_aras_natura2000_20190528:

path_natura2000 <- here::here(
  "data",
  "external",
  "protected_areas_natura2000_20190528"
)
if (!file.exists(path_natura2000)) {
  temp <- tempfile()
  download.file(
    "http://ftp.eea.europa.eu/www/natura2000/Natura2000_end2018_Shapefile.zip",
    temp
  )
  unzip(temp, exdir = path_natura2000)
  unlink(temp)
}

We read the shapefile:

protected_areas <- st_read(path_natura2000)

## Reading layer `Natura2000_end2018_epsg3035' from data source `C:\Users\damiano_oldoni\Documents\INBO\repositories\pipeline\data\external\protected_areas_natura2000_20190528' using driver `ESRI Shapefile'
## Simple feature collection with 27856 features and 5 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 746551.8 ymin: 905053.4 xmax: 6506130 ymax: 5298655
## CRS:            3035

protected_areas_crs <- st_crs(protected_areas)

Check whether grid and protected areas have same CRS:

be_grid_crs == protected_areas_crs

## [1] TRUE

Select Belgian protected areas

We are interested in protected areas in Belgium:

protected_areas <-
  protected_areas %>%
  filter(MS == "BE")

## filter: removed 27,546 rows (99%), 310 rows remaining

Summary by site type:

protected_areas %>%
  as_tibble() %>%
  group_by(SITETYPE) %>%
  count()

## group_by: one grouping variable (SITETYPE)

## count: now 3 rows and 2 columns, one group variable remaining (SITETYPE)

Read Belgian regions

The Belgian regions are available by package BelgiumMaps.StatBel:

data(BE_ADMIN_REGION)
BE_ADMIN_REGION <- st_as_sf(BE_ADMIN_REGION)

## Loading required package: sp

## Warning: package 'sp' was built under R version 3.6.3

Add regional information to Belgian proteted areas

In this section we add three columns to know whether a protected area belong to one or more Belgian regions.

Flemish protected areas

Add column flanders to indicate whether the protected area is a Flemish protected areas (TRUE/FALSE):

protected_areas_flanders <- tibble(
  SITECODE = c(
    "BE2200043",
    "BE2200032",
    "BE2100017",
    "BE2300006",
    "BE2500121",
    "BE2100024",
    "BE2400008",
    "BE2200030",
    "BE2200037",
    "BE2200041",
    "BE2200042",
    "BE2200038",
    "BE2101437",
    "BE2219312",
    "BE2500831",
    "BE2100323",
    "BE2100045",
    "BE2100020",
    "BE2200028",
    "BE2400009",
    "BE2200033",
    "BE2300044",
    "BE2400012",
    "BE2200626",
    "BE2200727",
    "BE2200029",
    "BE2400010",
    "BE2223316",
    "BE2301235",
    "BE2301336",
    "BE2300007",
    "BE2422315",
    "BE2524317",
    "BE2100026",
    "BE2100040",
    "BE2200031",
    "BE2500003",
    "BE2200034",
    "BE2100424",
    "BE2400014",
    "BE2101538",
    "BE2200039",
    "BE2500932",
    "BE2500001",
    "BE2218311",
    "BE2101639",
    "BE2217310",
    "BE2300222",
    "BE2400011",
    "BE2500002",
    "BE2220313",
    "BE2200035",
    "BE2100015",
    "BE2300005",
    "BE2301134",
    "BE2100016",
    "BE2100019",
    "BE2200036",
    "BE2500004",
    "BE2200525",
    "BE2501033",
    "BE2221314"
  )) %>%
  mutate(flanders = TRUE,
         wallonia = FALSE,
         brussels = FALSE)

## mutate: new variable 'flanders' with one unique value and 0% NA

##         new variable 'wallonia' with one unique value and 0% NA

##         new variable 'brussels' with one unique value and 0% NA

Brussels protected areas

Add column brussels to indicate whether the protected area is located in Brussels region (TRUE/FALSE):

protected_areas_brussels <- tibble(
  SITECODE = c("BE1000001",
               "BE1000002",
               "BE1000003")) %>%
  mutate(flanders = FALSE,
         wallonia = FALSE,
         brussels = TRUE)

## mutate: new variable 'flanders' with one unique value and 0% NA

##         new variable 'wallonia' with one unique value and 0% NA

##         new variable 'brussels' with one unique value and 0% NA

Protected areas under federal administration

Areas not belonging to any region as they are under federal administration:

protected_areas_federal <- tibble(
  SITECODE = c("BEMNZ0001",
               "BEMNZ0002",
               "BEMNZ0003",
               "BEMNZ0004",
               "BEMNZ0005")) %>%
  mutate(flanders = FALSE,
         wallonia = FALSE,
         brussels = FALSE)

## mutate: new variable 'flanders' with one unique value and 0% NA

##         new variable 'wallonia' with one unique value and 0% NA

##         new variable 'brussels' with one unique value and 0% NA

Wallonia protected areas

All the other protected areas are in Wallonia (wallonia = TRUE):

protected_areas_wallonia <- tibble(
  SITECODE = as.character(protected_areas$SITECODE[
    !protected_areas$SITECODE %in% c(protected_areas_flanders$SITECODE,
                                     protected_areas_brussels$SITECODE,
                                     protected_areas_federal$SITECODE)
  ])) %>%
  mutate(flanders = FALSE,
         wallonia = TRUE,
         brussels = FALSE)

## mutate: new variable 'flanders' with one unique value and 0% NA

##         new variable 'wallonia' with one unique value and 0% NA

##         new variable 'brussels' with one unique value and 0% NA

Add regional information

We can now merge the regional information

protected_areas_regional_info <-
  bind_rows(protected_areas_flanders,
            protected_areas_wallonia,
            protected_areas_brussels,
            protected_areas_federal) %>%
  mutate(SITECODE = as.factor(SITECODE))

## mutate: converted 'SITECODE' from character to factor (0 new NA)

to add it to protected_areas:

protected_areas <- 
  protected_areas %>%
  left_join(protected_areas_regional_info,
            by = "SITECODE")

## Warning: Column `SITECODE` joining factors with different levels, coercing to
## character vector

## left_join: added 3 columns (flanders, wallonia, brussels)

## Warning: Column `SITECODE` joining factors with different levels, coercing to
## character vector

## Warning: Column `SITECODE` joining factors with different levels, coercing to
## character vector

##            > rows only in x     0

##            > rows only in y  (  0)

##            > matched rows     310

##            >                 =====

##            > rows total       310

Preview:

protected_areas %>%
  as_tibble() %>%
  select(SITECODE, SITENAME, SITETYPE, flanders, wallonia, brussels) %>%
  head()

## select: dropped 3 variables (RELEASE_DA, MS, geometry)

Distribution protected areas

Number of protected areas per type:

protected_areas %>%
  st_drop_geometry() %>%
  group_by(SITETYPE) %>%
  count()

## group_by: one grouping variable (SITETYPE)

## count: now 3 rows and 2 columns, one group variable remaining (SITETYPE)

In Flanders:

protected_areas %>%
  filter(flanders == TRUE) %>%
  st_drop_geometry() %>%
  group_by(SITETYPE) %>%
  count()

## filter: removed 248 rows (80%), 62 rows remaining

## group_by: one grouping variable (SITETYPE)

## count: now 2 rows and 2 columns, one group variable remaining (SITETYPE)

# Transform to wgs84
protected_areas_wgs84 <- 
  protected_areas %>%
          st_transform(crs = 4326)
# Flemish protected areas
protected_areas_wgs84_flanders <- 
  protected_areas_wgs84 %>%
  filter(flanders == TRUE)

## filter: removed 248 rows (80%), 62 rows remaining

prot_area_fl_palette <- colorFactor(
  topo.colors(nrow(protected_areas_flanders)),
  protected_areas_flanders$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_flanders %>%
                filter(SITETYPE == "A"),
              fillColor = ~prot_area_fl_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_flanders %>%
                filter(SITETYPE == "B"),
              fillColor = ~prot_area_fl_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B"),
    options = layersControlOptions(collapsed = FALSE))

## filter: removed 38 rows (61%), 24 rows remaining

## filter: removed 24 rows (39%), 38 rows remaining

In Wallonia:

protected_areas %>%
  filter(wallonia == TRUE) %>%
  st_drop_geometry() %>%
  group_by(SITETYPE) %>%
  count()

## filter: removed 70 rows (23%), 240 rows remaining

## group_by: one grouping variable (SITETYPE)

## count: now 3 rows and 2 columns, one group variable remaining (SITETYPE)

protected_areas_wgs84_wallonia <- 
  protected_areas_wgs84 %>%
  filter(wallonia == TRUE)

## filter: removed 70 rows (23%), 240 rows remaining

prot_area_wallonia_palette <- colorFactor(
  topo.colors(nrow(protected_areas_wallonia)),
  protected_areas_wallonia$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                filter(SITETYPE == "A"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                filter(SITETYPE == "B"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addPolygons(data = protected_areas_wgs84_wallonia %>%
                filter(SITETYPE == "C"),
              fillColor = ~prot_area_wallonia_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type C") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B", "type C"),
    options = layersControlOptions(collapsed = FALSE))

## filter: removed 238 rows (99%), 2 rows remaining

## filter: removed 228 rows (95%), 12 rows remaining

## filter: removed 14 rows (6%), 226 rows remaining

In Brussels:

protected_areas %>%
  filter(brussels == TRUE) %>%
  st_drop_geometry() %>%
  group_by(SITETYPE) %>%
  count()

## filter: removed 307 rows (99%), 3 rows remaining

## group_by: one grouping variable (SITETYPE)

## count: now one row and 2 columns, one group variable remaining (SITETYPE)

protected_areas_wgs84_brussels <- 
  protected_areas_wgs84 %>%
  filter(brussels == TRUE)

## filter: removed 307 rows (99%), 3 rows remaining

prot_area_brussels_palette <- colorFactor(
  topo.colors(nrow(protected_areas_brussels)),
  protected_areas_brussels$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_brussels %>%
                filter(SITETYPE == "B"),
              fillColor = ~prot_area_brussels_palette(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type B"),
    options = layersControlOptions(collapsed = FALSE))

## filter: no rows removed

Areas under federal administration

protected_areas %>%
  filter(flanders == FALSE,
         wallonia == FALSE,
         brussels == FALSE) %>%
  st_drop_geometry() %>%
  group_by(SITETYPE) %>%
  count()

## filter: removed 305 rows (98%), 5 rows remaining

## group_by: one grouping variable (SITETYPE)

## count: now 2 rows and 2 columns, one group variable remaining (SITETYPE)

# Protected areas under federal administration
protected_areas_wgs84_federal <- 
  protected_areas_wgs84 %>%
  filter(flanders == FALSE,
         wallonia == FALSE,
         brussels == FALSE)

## filter: removed 305 rows (98%), 5 rows remaining

prot_area_palette_federal <- colorFactor(
  topo.colors(nrow(protected_areas_federal)),
  protected_areas_federal$SITECODE)
leaflet() %>%
  addTiles() %>%
  addPolygons(data = protected_areas_wgs84_federal %>%
                filter(SITETYPE == "A"),
              fillColor = ~prot_area_palette_federal(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type A") %>%
  addPolygons(data = protected_areas_wgs84_federal %>%
                filter(SITETYPE == "B"),
              fillColor = ~prot_area_palette_federal(SITECODE),
              fillOpacity = 0.7,
              color = "black",
              weight = 0.5,
              opacity = 1,
              label = ~SITENAME,
              popup = ~SITECODE,
              group = "type B") %>%
  addLayersControl(
    overlayGroups = c("type A", "type B"),
    options = layersControlOptions(collapsed = FALSE))

## filter: removed 2 rows (40%), 3 rows remaining

## filter: removed 3 rows (60%), 2 rows remaining

Define Bird and Habitat directive areas

Bird directive areas

Bird directive areas (SPA areas) are defined as the areas with SITETYPE: A and C.

protected_areas_bird <-
  protected_areas %>%
  filter(SITETYPE %in% c("A", "C"))

## filter: removed 55 rows (18%), 255 rows remaining

Habitat directive areas

Habitat directive areas (SAC areas) are defined as the protected areas with SITETYPE values B and C.

protected_areas_habitat <-
  protected_areas %>%
  filter(SITETYPE %in% c("B", "C"))

## filter: removed 29 rows (9%), 281 rows remaining

Intersect grid of Belgium with protected areas

Interset grid with protected areas

Intersect cells of the grid and protected areas:

grid_intersects_areas <- st_intersects(be_grid, protected_areas,
  sparse = FALSE
)
grid_intersects_areas <- as_tibble(grid_intersects_areas,
  .name_repair = "universal"
)

## New names:
## * `` -> ...1
## * `` -> ...2
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * ...

# Assign column names based on protected areas code
prot_areas_sitecode <- as.character(protected_areas$SITECODE)
names(grid_intersects_areas) <- prot_areas_sitecode
# Add 1x1km grid cellcode as new column
grid_intersects_areas <-
  grid_intersects_areas %>%
  mutate(CELLCODE = as.character(droplevels(be_grid$CELLCODE))) %>%
  select(CELLCODE, everything())

## mutate: new variable 'CELLCODE' with 51,726 unique values and 0% NA

Natura2000 protected areas

Select cells which intersect any of the protected areas:

cells_protected <-
  grid_intersects_areas %>%
  filter_at(vars(starts_with("BE")), any_vars(. == TRUE))

## filter_at: removed 38,988 rows (75%), 12,738 rows remaining

Get grid cells for each protected area

For each protected area we can get the list of intersecting cells:

cells_of_prot_areas <-
  map_dfr(
    prot_areas_sitecode,
    function(pa_cell_code) {
      site_type <-
        protected_areas %>%
        dplyr::filter(SITECODE == pa_cell_code) %>%
        dplyr::pull(SITETYPE)
      cells <-
        cells_protected %>%
        # avoid long logs by using dplyr functions instead of tidylog
        dplyr::select(CELLCODE, !!sym(pa_cell_code)) %>%
        dplyr::filter(!!sym(pa_cell_code) == TRUE) %>%
        dplyr::mutate(
          SITECODE = pa_cell_code,
          SITETYPE = as.character(droplevels(site_type))
        ) %>%
        dplyr::select(SITECODE, SITETYPE, CELLCODE)
    }
  )

For example, grid cells intersecting protected area BE1000003:

cells_of_prot_areas %>%
  filter(SITECODE == "BE1000003")

## filter: removed 14,806 rows (>99%), 8 rows remaining

Distribution of the number of cells covering the protected areas:

cells_of_prot_areas %>%
  group_by(SITECODE) %>%
  count() %>%
  arrange(desc(n)) %>%
  ggplot() +
  geom_histogram(aes(x = n), bins = 30) +
  scale_x_log10() +
  xlab("Number of cells (km2)") +
  ylab("Number of protected areas")

## group_by: one grouping variable (SITECODE)

## count: now 310 rows and 2 columns, one group variable remaining (SITECODE)

Natura2000

We can now define whether a cell of the reference grid interesects a protected area as defined by Natura2000. We add this information in column natura2000:

cellcode_protected_natura2000 <-
  cells_protected %>%
  as_tibble() %>%
  pull(CELLCODE)

be_grid_membership_protected_areas <-
  be_grid %>%
  mutate(natura2000 = ifelse(CELLCODE %in% cellcode_protected_natura2000,
    TRUE,
    FALSE
  ))

## mutate: new variable 'natura2000' with 2 unique values and 0% NA

Bird directive areas

We add also a column, spa, assessing whether the cell intersects a SPA area:

cells_protected_spa <-
  cells_protected %>%
  select(
    CELLCODE,
    one_of(as.character(protected_areas_bird$SITECODE))
  ) %>%
  filter_at(vars(starts_with("BE")), any_vars(. == TRUE)) %>%
  as_tibble() %>%
  pull(CELLCODE)

## select: dropped 55 variables (BE2200043, BE2100017, BE2200032, BE2100024, BE2400008, …)

## filter_at: removed 3,387 rows (27%), 9,351 rows remaining

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  mutate(spa = ifelse(CELLCODE %in% cells_protected_spa,
    TRUE,
    FALSE
  ))

## mutate: new variable 'spa' with 2 unique values and 0% NA

Habitat directive areas

We add also a column, habitat, assessing whether the cell intersects a habitat directive area:

cells_protected_hab <-
  cells_protected %>%
  select(
    CELLCODE,
    one_of(as.character(protected_areas_habitat$SITECODE))
  ) %>%
  filter_at(vars(starts_with("BE")), any_vars(. == TRUE)) %>%
  as_tibble() %>%
  pull(CELLCODE)

## select: dropped 29 variables (BE2500121, BE2101437, BE2219312, BE2500831, BE2100323, …)

## filter_at: removed 904 rows (7%), 11,834 rows remaining

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  mutate(habitat = ifelse(CELLCODE %in% cells_protected_hab,
    TRUE,
    FALSE
  ))

## mutate: new variable 'habitat' with 2 unique values and 0% NA

Summary

How many cells contain/intersect protected areas?

be_grid_membership_protected_areas <-
  be_grid_membership_protected_areas %>%
  as_tibble()
be_grid_membership_protected_areas %>%
  group_by(natura2000) %>%
  count()

## group_by: one grouping variable (natura2000)

## count: now 2 rows and 2 columns, one group variable remaining (natura2000)

How many cells contain/intersect SPA areas?

be_grid_membership_protected_areas %>%
  group_by(spa) %>%
  count()

## group_by: one grouping variable (spa)

## count: now 2 rows and 2 columns, one group variable remaining (spa)

How many cells contain/intersect habitat directive areas?

be_grid_membership_protected_areas %>%
  group_by(habitat) %>%
  count()

## group_by: one grouping variable (habitat)

## count: now 2 rows and 2 columns, one group variable remaining (habitat)

Save data

We save be_grid_membership_protected_areas as tab-separated file (tsv):

write_tsv(be_grid_membership_protected_areas,
  path = here::here(
    "data",
    "interim",
    "intersect_EEA_ref_grid_protected_areas.tsv"
  ),
  na = ""
)

We also save cells_of_prot_areas in the same format:

write_tsv(cells_of_prot_areas,
  path = here::here(
    "data",
    "interim",
    "EEA_ref_grid_cells_covering_protected_areas.tsv"
  ),
  na = ""
)

We are also interested in saving basic information about protected areas:

protected_areas_metadata <-
  protected_areas %>%
  st_drop_geometry() %>%
  select(SITECODE,
         SITENAME,
         SITETYPE,
         flanders,
         wallonia,
         brussels)

## select: dropped 2 variables (RELEASE_DA, MS)

protected_areas_metadata

We save these metadata:

protected_areas_metadata %>%
  write_tsv(
    path = here::here(
      "data",
      "interim",
      "protected_areas_metadata.tsv"
    ),
    na = ""
  )